https://jakubnowosad.com/posts/2024-10-20-spatcomp-bp2/

library(terra)
terra 1.7.23
ndvi2018_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2018_tartu.tif")
ndvi2023_tartu = rast("https://github.com/Nowosad/comparing-spatial-patterns-2024/raw/refs/heads/main/data/ndvi2023_tartu.tif")
plot(ndvi2018_tartu, main = "Tartu (2000)")

plot(ndvi2023_tartu, main = "Tartu (2018)")

ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
terra::plot(ndvi_diff)

plot_div = function(r, ...){
  r_range = range(values(r), na.rm = TRUE, finite = TRUE)
  max_abs = max(abs(r_range))
  new_range = c(-max_abs, max_abs)
  plot(r, col = hcl.colors(100, palette = "prgn"), range = new_range, ...)
}
plot_div(ndvi_diff)

ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
hist(ndvi_diff)

library(yardstick)
ndvi_rmse = rmse_vec(values(ndvi2023_tartu)[, 1], values(ndvi2018_tartu)[, 1])
ndvi_rmse
[1] 0.2191853
library(diffeR)
Loading required package: ggplot2
ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu)
[1] "The mean of grid1 is less than the mean of grid2"
ndvi_mad$Total
[1] 0.184306
ndvi_cor = focalPairs(c(ndvi2023_tartu, ndvi2018_tartu), w = 5,
                      fun = "pearson", na.rm = TRUE)
plot_div(ndvi_cor)

library(geodiv)

Attaching package: ‘geodiv’

The following objects are masked from ‘package:terra’:

    rotate, sds
window = matrix(1, nrow = 5, ncol = 5)
ndvi2018_tartu_sa_mw = focal_metrics(ndvi2018_tartu, window = window,
                               metric = "sa", progress = FALSE)
ndvi2023_tartu_sa_mw = focal_metrics(ndvi2023_tartu, window = window,
                               metric = "sa", progress = FALSE)
ndvi_diff_sa_mw = ndvi2023_tartu_sa_mw$sa - ndvi2018_tartu_sa_mw$sa
plot_div(ndvi_diff_sa_mw)

plot(ndvi2023_tartu)

plot(ndvi2023_tartu_sa_mw$sa)

library(rasterdiv)
ndvi2018_tartu_int = as.int(ndvi2018_tartu * 100)
ndvi2023_tartu_int = as.int(ndvi2023_tartu * 100)
ndvi2018_tartu_rao = paRao(ndvi2018_tartu_int, window = 5, progBar = FALSE)
Warning: Simplify=0. Rounding data to 0 decimal places.

Processing alpha: 1 Moving Window: 5


Processing alpha: 1 Moving Window: 5
Error in utils::combn(p, m = 2, FUN = prod, na.rm = TRUE) : n < m
ndvi_diff = ndvi2023_tartu - ndvi2018_tartu
ndvi_diff_autocor = autocor(ndvi_diff, method = "moran", global = FALSE)
plot_div(ndvi_diff, main = "Difference")

plot_div(ndvi_diff_autocor, main = "Moran's I of the difference")

library(SSIMmap)
ndvi_ssim = ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = FALSE, w = 3)
plot_div(ndvi_ssim[[1]])

library(waywiser)
cell_sizes = seq(50, 300, by = 50)
ndvi_multi_scale = ww_multi_scale(truth = ndvi2018_tartu, estimate = ndvi2023_tartu,
                                 metrics = list(yardstick::rmse), 
                                 cellsize = cell_sizes,
                                 progress = FALSE)
ℹ The package "exactextractr" is required.
✖ Would you like to install it?

1: Yes
2: No
1
Installing package into ‘C:/Users/stahlm/AppData/Local/R/win-library/4.2’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/bin/windows/contrib/4.2/exactextractr_0.10.0.zip'
Content type 'application/zip' length 1448767 bytes (1.4 MB)
downloaded 1.4 MB
package ‘exactextractr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\stahlm\AppData\Local\Temp\RtmpWCibE9\downloaded_packages
Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.Warning: No CRS specified for polygons; assuming they have the same CRS as the raster.
ndvi_multi_scale
library(diffeR)
ndvi_mad = MAD(ndvi2023_tartu, ndvi2018_tartu, eval = "multiple")
[1] "The mean of grid1 is less than the mean of grid2"

ndvi_mad
library(SSIMmap)
ssim_raster(ndvi2018_tartu, ndvi2023_tartu, global = TRUE)
SSIM: 0.63845 SIM: 0.90432 SIV: 0.90778 SIP: 0.75671 
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpodHRwczovL2pha3Vibm93b3NhZC5jb20vcG9zdHMvMjAyNC0xMC0yMC1zcGF0Y29tcC1icDIvDQoNCmBgYHtyfQ0KbGlicmFyeSh0ZXJyYSkNCm5kdmkyMDE4X3RhcnR1ID0gcmFzdCgiaHR0cHM6Ly9naXRodWIuY29tL05vd29zYWQvY29tcGFyaW5nLXNwYXRpYWwtcGF0dGVybnMtMjAyNC9yYXcvcmVmcy9oZWFkcy9tYWluL2RhdGEvbmR2aTIwMThfdGFydHUudGlmIikNCm5kdmkyMDIzX3RhcnR1ID0gcmFzdCgiaHR0cHM6Ly9naXRodWIuY29tL05vd29zYWQvY29tcGFyaW5nLXNwYXRpYWwtcGF0dGVybnMtMjAyNC9yYXcvcmVmcy9oZWFkcy9tYWluL2RhdGEvbmR2aTIwMjNfdGFydHUudGlmIikNCnBsb3QobmR2aTIwMThfdGFydHUsIG1haW4gPSAiVGFydHUgKDIwMDApIikNCnBsb3QobmR2aTIwMjNfdGFydHUsIG1haW4gPSAiVGFydHUgKDIwMTgpIikNCmBgYA0KDQpgYGB7cn0NCm5kdmlfZGlmZiA9IG5kdmkyMDIzX3RhcnR1IC0gbmR2aTIwMThfdGFydHUNCnBsb3QobmR2aV9kaWZmKQ0KYGBgDQoNCmBgYHtyfQ0KcGxvdF9kaXYgPSBmdW5jdGlvbihyLCAuLi4pew0KICByX3JhbmdlID0gcmFuZ2UodmFsdWVzKHIpLCBuYS5ybSA9IFRSVUUsIGZpbml0ZSA9IFRSVUUpDQogIG1heF9hYnMgPSBtYXgoYWJzKHJfcmFuZ2UpKQ0KICBuZXdfcmFuZ2UgPSBjKC1tYXhfYWJzLCBtYXhfYWJzKQ0KICBwbG90KHIsIGNvbCA9IGhjbC5jb2xvcnMoMTAwLCBwYWxldHRlID0gInByZ24iKSwgcmFuZ2UgPSBuZXdfcmFuZ2UsIC4uLikNCn0NCnBsb3RfZGl2KG5kdmlfZGlmZikNCmBgYA0KDQpgYGB7cn0NCm5kdmlfZGlmZiA9IG5kdmkyMDIzX3RhcnR1IC0gbmR2aTIwMThfdGFydHUNCmhpc3QobmR2aV9kaWZmKQ0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeSh5YXJkc3RpY2spDQpuZHZpX3Jtc2UgPSBybXNlX3ZlYyh2YWx1ZXMobmR2aTIwMjNfdGFydHUpWywgMV0sIHZhbHVlcyhuZHZpMjAxOF90YXJ0dSlbLCAxXSkNCm5kdmlfcm1zZQ0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeShkaWZmZVIpDQpuZHZpX21hZCA9IE1BRChuZHZpMjAyM190YXJ0dSwgbmR2aTIwMThfdGFydHUpDQpgYGANCmBgYHtyfQ0KbmR2aV9tYWQkVG90YWwNCmBgYA0KDQpgYGB7cn0NCm5kdmlfY29yID0gZm9jYWxQYWlycyhjKG5kdmkyMDIzX3RhcnR1LCBuZHZpMjAxOF90YXJ0dSksIHcgPSA1LA0KICAgICAgICAgICAgICAgICAgICAgIGZ1biA9ICJwZWFyc29uIiwgbmEucm0gPSBUUlVFKQ0KcGxvdF9kaXYobmR2aV9jb3IpDQpgYGANCg0KYGBge3J9DQpsaWJyYXJ5KGdlb2RpdikNCndpbmRvdyA9IG1hdHJpeCgxLCBucm93ID0gNSwgbmNvbCA9IDUpDQpuZHZpMjAxOF90YXJ0dV9zYV9tdyA9IGZvY2FsX21ldHJpY3MobmR2aTIwMThfdGFydHUsIHdpbmRvdyA9IHdpbmRvdywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRyaWMgPSAic2EiLCBwcm9ncmVzcyA9IEZBTFNFKQ0KbmR2aTIwMjNfdGFydHVfc2FfbXcgPSBmb2NhbF9tZXRyaWNzKG5kdmkyMDIzX3RhcnR1LCB3aW5kb3cgPSB3aW5kb3csDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbWV0cmljID0gInNhIiwgcHJvZ3Jlc3MgPSBGQUxTRSkNCm5kdmlfZGlmZl9zYV9tdyA9IG5kdmkyMDIzX3RhcnR1X3NhX213JHNhIC0gbmR2aTIwMThfdGFydHVfc2FfbXckc2ENCnBsb3RfZGl2KG5kdmlfZGlmZl9zYV9tdykNCmBgYA0KDQoNCmBgYHtyfQ0KcGxvdChuZHZpMjAyM190YXJ0dSkNCnBsb3QobmR2aTIwMjNfdGFydHVfc2FfbXckc2EpDQpgYGANCg0KYGBge3J9DQpsaWJyYXJ5KHJhc3RlcmRpdikNCm5kdmkyMDE4X3RhcnR1X2ludCA9IGFzLmludChuZHZpMjAxOF90YXJ0dSAqIDEwMCkNCm5kdmkyMDIzX3RhcnR1X2ludCA9IGFzLmludChuZHZpMjAyM190YXJ0dSAqIDEwMCkNCm5kdmkyMDE4X3RhcnR1X3JhbyA9IHBhUmFvKG5kdmkyMDE4X3RhcnR1X2ludCwgd2luZG93ID0gNSwgcHJvZ0JhciA9IEZBTFNFKQ0KbmR2aTIwMjNfdGFydHVfcmFvID0gcGFSYW8obmR2aTIwMjNfdGFydHVfaW50LCB3aW5kb3cgPSA1LCBwcm9nQmFyID0gRkFMU0UpDQpuZHZpX3Jhb19kaWZmID0gbmR2aTIwMjNfdGFydHVfcmFvW1sxXV1bWzFdXSAtIG5kdmkyMDE4X3RhcnR1X3Jhb1tbMV1dW1sxXV0NCnBsb3RfZGl2KG5kdmlfcmFvX2RpZmYpDQpgYGANCg0KYGBge3J9DQpuZHZpX2RpZmYgPSBuZHZpMjAyM190YXJ0dSAtIG5kdmkyMDE4X3RhcnR1DQpuZHZpX2RpZmZfYXV0b2NvciA9IGF1dG9jb3IobmR2aV9kaWZmLCBtZXRob2QgPSAibW9yYW4iLCBnbG9iYWwgPSBGQUxTRSkNCnBsb3RfZGl2KG5kdmlfZGlmZiwgbWFpbiA9ICJEaWZmZXJlbmNlIikNCnBsb3RfZGl2KG5kdmlfZGlmZl9hdXRvY29yLCBtYWluID0gIk1vcmFuJ3MgSSBvZiB0aGUgZGlmZmVyZW5jZSIpDQpgYGANCg0KDQoNCmBgYHtyfQ0KbGlicmFyeShTU0lNbWFwKQ0KbmR2aV9zc2ltID0gc3NpbV9yYXN0ZXIobmR2aTIwMThfdGFydHUsIG5kdmkyMDIzX3RhcnR1LCBnbG9iYWwgPSBGQUxTRSwgdyA9IDMpDQpwbG90X2RpdihuZHZpX3NzaW1bWzFdXSkNCmBgYA0KYGBge3J9DQpsaWJyYXJ5KHdheXdpc2VyKQ0KY2VsbF9zaXplcyA9IHNlcSg1MCwgMzAwLCBieSA9IDUwKQ0KbmR2aV9tdWx0aV9zY2FsZSA9IHd3X211bHRpX3NjYWxlKHRydXRoID0gbmR2aTIwMThfdGFydHUsIGVzdGltYXRlID0gbmR2aTIwMjNfdGFydHUsDQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBtZXRyaWNzID0gbGlzdCh5YXJkc3RpY2s6OnJtc2UpLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGNlbGxzaXplID0gY2VsbF9zaXplcywNCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHByb2dyZXNzID0gRkFMU0UpDQpuZHZpX211bHRpX3NjYWxlDQpgYGANCg0KYGBge3J9DQpsaWJyYXJ5KGRpZmZlUikNCm5kdmlfbWFkID0gTUFEKG5kdmkyMDIzX3RhcnR1LCBuZHZpMjAxOF90YXJ0dSwgZXZhbCA9ICJtdWx0aXBsZSIpDQpgYGANCmBgYHtyfQ0KbmR2aV9tYWQNCmBgYA0KDQoNCg0KYGBge3J9DQpsaWJyYXJ5KFNTSU1tYXApDQpzc2ltX3Jhc3RlcihuZHZpMjAxOF90YXJ0dSwgbmR2aTIwMjNfdGFydHUsIGdsb2JhbCA9IFRSVUUpDQpgYGANCg0K